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SECTION ONE 


Introduction 


1.0 OVERVIEW 

The design and analysis of even a simple multirate compensator can be a complex task. In this document we 
describe some Matlab M -Files which aid in multirate design. They include M-Files for modeling multirate systems, 
for computing optimum values of a multrrate compensator’s gains, and for generating time domain simulations. 

The remainder of the document is divided into three sections. In Section Two we review the basics of multirate 
modeling and our optimization algorithm. We also present the key concepts and notation which are utilized in 
Section Three. Section Three describes each M-File, detailing its inputs and outputs. The M-Files in this section 
are divided into three categories: modeling, simulation, and synthesis. Finally, in Section Four we conclude with a 
multirate design example. 

1.1 SOFTWARE REQUIREMENTS 

In addition to the standard Matlab toolbox routines, this software uses M-Files from both the Control_Toolbox 
and the Optimization_Toolbox. These two toolboxes must be present for the multirate software to operate. 
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SECTION TWO 

Background 


2.0 OVERVIEW 

In the following paragraphs we present the notation used in the remaining sections and review some key 
concepts which will be helpful to those using this software. A detailed explanation of the theory behind our M-Files 
can be found in References 1-5. 

2.1 A MULTIRATE FEEDBACK SYSTEM 

A multirate sampled-data system consists of a continuous plant in feedback with a multirate compensator. A 
block diagram of such a system is shown in Figure 2.1. The vector y s in this figure represents the continuous plant 

sensor output. Each element of y s can be sampled at an independent rate. The vector y represents the sampled 
value of y s available to the digital processor. (In our double index notation, the discrete signal p(m y n) results from 
sampling the continuous signal p(t) at the times t = ( mN + n)T\ where the integer N is the period of repetition; T is 
the sampling time; m =0, 1, ... ; and n- 0, 1, ..., AM.) The digital processor obtains the current value of y and 
combines it with the current processor state vector, z, using the state space structure shown in Figure 2.1. Each 
element of the processor state vector, z , can be updated at an independent rate. The continuous output from the 
compensator, represented by the vector u , is formed by holding the output from the digital processor, w, with a zero- 
order-hold. Each element of the vector u can be held at an independent rate to form u. The vectors v and w 
represent the discrete sensor and continuous process noise respectively. 

Conceptually, one can divide the multirate compensator into two parts, the “sampling schedule” and the digital 
processor gains. This is the approach used in the synthesis and analysis software. The “sampling schedule” is a 
description of when each compensator input is sampled and when each output and processor state is updated, while 
the digital processor gains determine the dynamics of the digital processor. In the following paragraphs we discuss 
each in detail. 



Figure 2.1 Closed-loop Multirate System 
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2.1.1 The Sampling Schedule 

In general, the sampling and updating of y Si z, and u in Figure 2. 1 can occur at any time. We, however, require 
that these sample and update activities occur only at integer multiples of some fixed time, called the shortest time 
period (STP). The actual value of the S TP is arbitrary, but is often a function of the hardware and software used to 
implement the control law. We also require that the sampling and updating activities of the sensors, states and 
outputs repeat themselves after some fixed period of time. The period of repetition of the sampling schedule is 
called the basic time period (BTP). Finally, we define 

"DTp 

the integer N - r=— and the value T = STP (2. 1 ) 

olr 

In our double index notation, the first index (m) in p(m,n) indicates the integer number of BTP’s which have elapsed 
when the sample/update occurred and the second index (n) indicates the integer number of STP’s which have 
elapsed within the current BTP when the sample/update occurred. 

We can represent the sampling schedule for the multirate compensator graphically, as shown in Figure 2.2. 
The figure shows a time line for each sampler, processor state, and zero-order-hold. The time line is divided into 
one STP increments. On the left side of the time line is a description of the signal or state, being sampled or 
updated. On the right side is a description of the particular activity represented by the time line, e.g. state update, 
sampler, or zero-order-hold. Circles on each time line indicate when a sample or update activity associated with that 
particular signal or state takes place. Usually the sampling schedule is shown for only one BTP since the sampling 
schedule repeats itself every BTP. 

In most applications, the sampling/updating activities for a given sensor, output or state will be periodic within 
the BTP, as is shown in Figure 2.2. However, the sampling/updating activities do not have to be periodic within the 
BTP. The only requirement is that the samphng/updating activities have some period of repetition (the BTP) and 
that they occur at integer multiples of the STP. Once the STP and BTP have been selected, the designer can 
arbitrarily specify sampling/updating activities at any multiple of the STP with in one BTP. An example of a 
multirate sampling schedule in which the sampling/updating activities are not periodic within the BTP is shown in 
Figure 2.3. A sampling policy like this might be used to multiplex multiple inputs through a single analog to digital 
converter. 

2.1.2 Digital Processor Gains 

The processor gains are the values of the matrices A 2 , B z , C 2 , and D z in Figure 2.1. Like the sampling schedule, 
they can be periodically time-varying with a period of repetition of one BTP. Generally, these matrices are free 
design parameters which can be adjusted by the designer to improve the performance of the multirate compensator. 
The synthesis software discussed later in this document can be used to calculate optimum values for these gains. 

22 THE EQUIVALENT TIME -INVARIANT SYSTEM (ETIS) 

A multirate compensator with the structure discussed in Section 2.1 can be modeled as a periodically time- 
varying single-rate compensator by appending appropriate hold states to the digital processor model. This new 
compensator has the form 

jt(m,7H-l) = A(n)x(m,n) + B(n)u(m,ri) (2.2a) 

y s (m y n) - C(n)x(m,ri) + D(n)u(mji) (2.2b) 

for n - 0,1,. .JV-1, and m - 0,1,... 

This compensator has a sampling period of one STP and a period of repetition of one BTP (or NT), y s (tn,n) 
represents the values of y(t ) sampled every STP; u(t) is formed by holding u(m t n) with a zero- order-hold. 
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Figure 2.2. Example multirate sampling schedule with periodic sampling/updating activity 
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Figure 2.3. Example multirate sampling schedule with aperiodic sampling activity 


The periodically time-varying single-rate compensator can be further transformed into a single-rate Equivalent 
Time-Invariant System (ETIS) with the form shown below. 

jt(rn+l,0) = Agx(m y 0) + B gup{m,0) (2.3a) 

>>£(m,0) = C£t(m,0) + DpUjAmfi) (2.3b) 

where 


and yp(m, 0) 




u(m ,0) 

y s (m> 1) 

utfm, 0) = 

u(m, 1) 

-y s (.mff- 1)_J 


-w(m //-!)- 


(2.4) 


We use the subscript E to denote vectors and matrices strictly associated with the ETIS. Notice that the ETTS 
input/output vectors are composite vectors containing the input/output values of the multirate (or periodically-time 
varying) system at N sampling times. Consequently, an ETIS is always MIMO even if the original system is SISO. 
If the multirate system has p inputs, q outputs and a sampling period of one STP then the ETTS is a single-rate linear 
time-invariant system with Np inputs, Nq outputs and a sampling period of one BTP. 

The ETIS is fundamental to the analysis of multirate systems. It allows one to evaluate the performance and 
stability of complex systems comprised of multirate, periodically time- varying and/or single-rate components using 
only modified linear time-invariant single-rate techniques. For example, to evaluate the stability of the system in 
Figure 2.1, we would first transform the multirate compensator into its ETIS with a given value for N. Then we 
would discretize the plant at the STP of the compensator using a zero- order-hold and transform the resulting single- 
rate system into an ETIS using the BTP of the compensator. Next, the plant and compensator ETIS s could be 
combined in feedback just as if they were traditional single-rate systems. Finally, we could determine the stability 
of the original multirate sampled-data system from the eigenvalues of its closed-loop ETIS. 
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23 OPTIMIZING THE DIGITAL PROCESSOR GAINS 

The synthesis software discussed later calculates optimum values of the digital processor gains, A z , B z , C z , and 
Z) z , by minimizing a quadratic cost function which reflects the performance of the closed-loop system in Figure 2.1 . 
The cost function has the form 


7 = lim E< 

/— > oo 


f T,.s 

T/.xl 

ra 

M' 

>c«]l 

[Jfc M 

U (0J 

1 m t 

Qi_ 



(2.5) 


where J is the cost associated with the closed-loop system shown in Figure 2.1. The vector y c is the continuous 
criterion output and w is the continuous control input Q \ , £>2 and M are the cost function weighting matrices which 
are selected by the designer. 

This cost function has the same form as that minimized by a continuous time LQR design. Thus the cost 
associated with the optimized multirate compensator and that of an LQR design can be compared directly. The 
designer can also use this fact to help select appropriate values for Q \ , 02 and AL 

To improve the robustness of the compensator, the optimization algorithm can optimize the digital processor 
gains for several different plant perturbations simultaneously. The resulting compensator will stabilize the plant at 
each perturbation and provide overall optimum performance. This is accomplished by minimizing the new cost 
function of Eqn. (2.6) which is the sum of the costs associated with each plant perturbation. 


J = 





(2.6) 


Here J J is the cost associated with the i ^ plant perturbation and there are Np plant perturbations. 
Optimum values of A z , B z , C z , and Z) z , occur when 


dJ 


= o, 



7 ^-°’ and izj: =0 


(2.7) 


Our algorithm uses a gradient type numerical search to determine values of the digital processor gains such that the 
conditions in Eqn. (2.7) are satisfied. Because the synthesis software uses an iterative process to determine optimum 
values for the digital processor gains, the user must provide the software with an initial guess for A z , B z> C z , and D z . 
The compensator corresponding to these values must stabilize every plant perturbation considered in (2.6). 
Obtaining a suitable stabilizing initial guess can sometimes be a difficult problem. We refer the interested reader to 
Reference L 
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Section Three 

M-File Documentation 


3.0 STANDARD VARIABLE DEFINITIONS 

Many of the M-Files discussed in this document require similar input variables or provide similar outputs. To 
simplify the documentation of these M-Files a set of standard variables are used throughout this document. They are 
defined in Table 3.1 with Matlab variables and functions bolded. 


Table 3.1 Standard Variable Definitions 


Variable 


Description 


pit 


nplt 

cmp 


Plant matrices in the form pit = [F, G; H, J] where the state-space representation of the 
plant is 



x(t) = Fx(t)+Gu(t) 

(3.1a) 


y(t) = Hx(t)+Ju(t) 

(3.1b) 

or 

x(m,n+l) = Fx(m,n ) + Gu{m,n) 

(3.1c) 


y(m,n) = Hx(m/t)+ Ju(m,n) 

(3. Id) 


depending on the value of stp defined later. 

Number of states in (3.1), or equivalently the number of rows in F 


Multirate compensator gain matrices. Given the state-space representation of the digital 
processor 

z(m,n + 1) = ^(n)z(m, n) + B z (n)y(m,n ) (3.2a) 

u(m,n) = C z (n)z(m,rt)+ D z (n)y{m y n) (3.2b) 


where m=0,l,... and n=0,l»—*N-l, so that the gain matrices, A z , B z , C z and Z) z , are 
periodically time- varying, cmp has the form 


cmp = 


T^(0) B z (oy 

'Ae(I) w 


~M N ~ 1) ^((V-DT 

[c z (0) D,(0)| 

c z ( 1) D z ( 1) 

9 y 

C z (N-l) Z) z (A^-l)J_ 


(3.3) 


If A z , F z , C z , and D z are constant, then 


A z (0) B z ( 0) 
C z ( 0) D z { 0) 


(3.4) 


and it is assumed that A Z (0)=A Z (1)= ... =A z (A r -l), F z (0)=2? z (l)= ... =B Z (N-1), etc. 

The software automatically deduces from the size of cmp, the value of ncmp and the 
number of compensator inputs whether the digital processor gains are periodically time- 
varying or not. 


Continued on following page... 
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Table 3.1 Standard Variable Definitions ( continued from previous page) 


Variable 


Description 


ncmp Number of digital processor states in Eqn. (3.2), or equivalently the number of rows in A z . 


stp 


stppbtp 

su 


If stp > 0 then stp is STP, the shortest sample/update period of the compensator. See 
Section 2.1.1, and pit describes a continuous plant. 

If stp = 0 then the plant matrices pit describe a discrete system whose sampling rate is 
one STP of the compensator, and pit describes a discrete plant. 

Number of STP’s per BTP. stppbtp = N. See Section 2.1.1. 

Sampling schedule for the compensator output, su has two forms 

Case I: su is a matrix, su must have as many columns as there are compensator outputs, 
and must have N rows (N = BTP/STP). Each element of su has a value of 1 or 0. A 1 
in the i^ 1 column and j^ 1 row of su indicates that the ft 1 compensator output is updated 
at the f 1 STP within the current BTP. A 0 indicates no update and the previous value 
is held. For example, given a two output system with the following sampling schedule 


Output #1 % — | — ? ZOH 
Output #2 — £ — ^ ZOH 

stp — J L_L 

^Z\ — BTP 


? = update activity 
A/= 3 


the corresponding value of su is 


su = 


1 0 | 
0 1 
1 1 


(3.5) 


Case II: su is a row vector. The i* column of su specifies the sampling period of the ft 
compensator output hold in number of STP’s. For example, su=[l 4] specifies that the 
first output is updated every STP and the second output is updated every fourth STP. 
This form assumes that the updating is synchronized on the first STP of the BTP. 

sy Sampling schedule for the compensator input, sy has the same form as su except that it 

specifies when the plant output (the compensator input) is sampled. 

sz Sampling schedule for the compensator states, sz has the same form as su except that it 

specifies when the compensator’ s digital processor states are updated. 
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3.1 MODELING M-FILES 

3.1.1 mr2etis 

Format : [ae ,be ,ce ,de )= mr2etis( cmp, ncmp , su, sy , sz , stppbtp) 

Description: Converts a multirate system into an ETIS with the form: 




' y(m, 0) 

x(m + 1, 0) = a er(m, 0) + be< 

y(m,n + 1) 

; 



y(m,N-l) 

w(m,0) 


y(m, 0) 

w(m,n + 1) 

► = cex(m,0) + de< 

y(m,n+ 1) 

u(m,N-l ) 




Inputs : cmp, ncmp , su , sy , sz, stppbtp The compensator description. See Table 3.1. 

Outputs : ae, be , ce , de The ETTS matrices in Eqn. (3.6) 


3.1.2 mr2ptv 

Format : [P, nu , ny, nz]=mr2ptv (cmp, ncmp , su , sy , sz, stppbtp ) 

Description : Converts a multirate system into a periodically time- varying system of the form 

x(m,n +1) = A(n)x{m > ri) + B(ri)y(m y ri) 
u(m y n)= C(n)r(m,n) + D(n)y(m,h) 


Inputs : cmp, ncmp , su , sy , sz, stppbtp The compensator description. See Table 3.1. 

Outputs: P The periodically time- varying system matrices in Eqn. (3.7), where 


rrD(o) c(0)i 

’D(l) C(l)‘ 


'D(N-l) C(tf-1)T 

o 

o 

.fi(l) A( 1)_ 

’ * 

[5(^-1) A(AT-1)J_ 


nu Number of inputs u in Eqn. (3.7) 

ny Number of outputs y in Eqn. (3.7) 

nz Number of states x in Eqn. (3.7) 


(3.6a) 


(3.6b) 


(3.7a) 

(3.7a) 


(3.8) 
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3.1.3 ptv2etis 

Format : [ae , be , ce , de]=ptv2etis(P, nu, ny , nz , stppbtp ) 

Description : Converts a periodically time-vaiying system into its ETTS of the form of Eqn. (3.6) 

Inputs : P, nu, ny, nz The periodically time- varying system matrices. See Section 3.L2 

stppbtp See Table 3 . 1 

Outputs : ae, be, ce, de The ET1S matrices. See Eqn. (3.6) 


3.1.4 sr2etis 

Format: [ae, be, ce, de}=sr2etis(a, b, c, d, stppbtp) 

Description. Converts the single-rate system in Eqn. (3.9) into an ETIS with an N of stppbtp, and the form of 
Eqn. (3.6) 

x(m,n +1) = a(ji)x(m,ri) + b{n)y{m,n) (3.9a) 

w(m,n)= c(n)x(m,n) + d(n)y(m,n ) (3.9b) 

Inputs: a, b, c, d The single-rate systems matrices. See Eqn. (3.9) 

stppbtp The desired N of the new ETIS. See Table 3.1 

Outputs: ae, be,ce,de The ETTS matrices. See Eqn. (3.6) 


3.1.5 stack 

Format : [y>= stack (w, stppbtp ) 

Description: Converts a traditional discrete- time vector w into an ETTS vector y 

Inputs: w A matrix whose columns contain the values of the vector y(m,n) for Nk successive samples times, 

where k is an integer. The matrix w has the form 

W = [y(0,0),y(0,l), ... ,y(0,A-l)y(l,0), y(l,l), ..., tfl.AM), .... y(k-hN-l)] (3.10) 

y is written using the double index notation of Section 2.1. 
stppbtp See Table 3.1. 

Outputs : y An ETTS vector of the values of yp (the ETTS of y) with the form 




3.1.6 


unstack 


Format : 
Description: 
Inputs : 

Outputs : 


[w)=unstack(y, stppbtp) 

Converts an ETIS vector y into a traditional discrete-time vector w. 

y An ETIS vector of the form of Eqn. (3.11) 
stppbtp See Table 3.1 

w A vector of the values of y(m,n) as in Eqn. (3.10) 


NOTE: Additional M-Files used in support of those described in this section are described in Appendix B. 
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3.2 SIMULATION M-FILES 

3.2.1 mrsim 

Format : [y, t] =mrsim(plt, nplt, cmp, ncrap, su, sy, sz, stppbtp , simpstp, inputl, outputl, ur , stp, XO) 

Description : Generates a time response simulation of a continuous or discrete plant in feedback with a multirate 

compensator. The routine assumes positive feedback as shown below with the continuous inputs 
and outputs of interest sampled or updated every stp /simpstp seconds. This configuration allows 
the user to simulate the inter- sample response of the multirate system. 



T f = stp/simpstp 

Inputs : pit, nplt A description of the plant, see Table 3.1. 

cmp, ncmp , su , sy , sz, stppbtp ,A description of the multirate compensator, see Table 3.1. 

simpstp This routine can calculate the response of the closed-loop system between the 
sample/update activities when the plant is continuous (stp > 0). simpstp is the number of 
simulation steps per STP of the compensator, simpstp must be an integer greater than zero and 
simpstp = 1 if stp = 0. 

inputl A vector specifying which plant inputs are connected to the compensator outputs. For 
example inputl = [14] indicates that the compensator outputs are connected to the first and 
fourth plant inputs. 

outputl A vector specifying which plant outputs are connected to the compensator inputs. For 
example outputl = [2 4] indicates that the compensator inputs are connected to the second and 
fourth plant outputs. 

ur The discrete input vector containing the values of the input at every (stp /simpstp) time for 
stppbtp simpstp sample times, ur has the form 

ur = [up(0,0), u^(0,l), ... , Up (£-1, stppbtp simpstp - 1)] (3.12) 

Note that in Eqn. (3.12) the sample/update period is one stp/simpstp. This sample period should 
not be confused with the sample period of the multirate compensator which is simply one stp. 
The second index in the double index notation in Eqn. (3.12) takes on values between zero and 
simpstp- stppbtp. 

stp see Table 3.1. 

XO A vector containing the initial conditions of the plant state vector. If XO is omitted, the initial 
conditions are assumed to be zero. 


i 


Outputs : 


y A vector containing the response of the plant every (stp/simpstp) time, y has the form 
y = [>p(0,0), y p (0,l), ... ,y p (k-l, stppbtp simpstp - 1)] 


(3.13) 
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Note that the second index in the double index notation corresponds to the second index of 
Eqn. (3.12) 

t a vector containing the times corresponding to the columns of y. When stp = 0 then t indicates the 
number of sampling periods. 

Comments : This M-File uses the ETIS of the system to compute the time domain response. Hence, if k in 

Eqn. (3.12) is not a multiple of simpstp jV where (N=BTP/STP) the system response is only 
simulated to the nearest multiple of simpstp N 

3.2.2 mrfeedback 

Format : [ae, be, ce, de] =mrfeedback( pit, nplt, cmp, ncmp, su, sy, sz,stppbtp, inputl, outputl , stp) 

Description : Creates a closed-loop ETIS system from a discrete or continuous plant, depending on stp , and a 

multirate compensator using positive feedback. 

Inputs: pit, nplt The plant description. See Table 3.1. 

cmp, ncmp, su, sy, sz, stppbtp The multirate compensator description. See Table 3.1. 
inputl A vector specifying which plant inputs are connected to the compensator outputs. For 
example inputl =[14] indicates that the compensator outputs are connected to the first and 
fourth plant inputs. 

outputl A vector specifying which plant outputs are connected to the compensator inputs. For 
example outputl = [2 4] indicates that the compensator inputs are connected to the second and 
fourth plant outputs. 

stp See Table 3.1. If stp is omitted mrfeedback assumes stp=0. 

ae, be, ce, de The system matrices of the ETIS closed-loop system with the form of Eqn. (3.6). The 
inputs and outputs of this system are the inputs and outputs of the discrete plant 


Outputs : 
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33 SYNTHESIS M-FILES 

The synthesis algorithm is implemented as a series of script files. They allow the designer to optimize the 
performance of the closed-loop system shown in Figure 3.1 by calculating values of A z , B Z ,C Z , and D z , such that the 
cost function of Eqn. (2.6) is minimized. There are three distinct steps to the optimization process: 1) entering the 
data which describes the problem; 2) discretizing and preprocessing the data; and 3) optimizing the compensator 
gain values. Each step of the optimization process is described in detail in the following paragraphs. 

Although the synthesis algorithm was designed to solve the multirate sampled-data problem, it is fairly general 
and can be applied to a variety of related problems as well. It can compute optimum compensator gains for a 
sampled-data system consisting of a continuous plant in feedback with a discrete compensator as shown in 
Figure 2. 1, or for a discrete system consisting of a discrete plant in feedback with a discrete compensator. In either 
case, the compensators may be single-rate or multirate and may have either time-invariant or periodically time- 
varying digital processor gains. 

33.1 Input Variables 

The synthesis algorithm looks for specific Matlab “workspace” variables to define the optimization problem. 
(A “workspace” variable is defined from the » prompt as opposed to being passed as an argument to a function.) 
These variables define things such as the plant dynamics, the cost function and the compensator structure and are 
defined in Table 3.2. The user must assign values to these variables before beginning the optimization. This can 
either be done manually or automatically from a script. An example of a script which defines these variables is 
provided in Appendix D. 



Digital Processor v<m;r7) 


Figure 3.1 Closed-loop Multirate System 
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Table 3.2 Matlab Workspace Variables Defining Optimization Problem 


Description 

The number of plant perturbations. Np=l indicates only the nominal pant is defined. See 
Section 2.3 

State-space matrices of the plant(s). Given the ft 1 continuous plant perturbation 

*(0 = ^•(r) + G,j^ ) ) J (3.14a) 

I yd( 3 = o- ] > 4b ) 


(3.14c) 

(3.14d) 


or the discrete plant perturbation, depending on the value of stp defined later, 

X((m, n + 1) = F/Jt; (m, n) + G/u^m, n) + w, (m, n) 


y c i(m,n) 

y s i(m,n) 


= H i x i (m,n) + Ji 


fui(m,nj 

0 


where u is the control input; W( is zero mean white noise input.; y c ( is the criterion output; 
and y s j are the outputs which are sampled by the compensator. Then 

pltJ \ F 1 ^]..P Np ^ P l] (3.15) 

/H [H 2 ^2 J ^Np ^Np 


For convenience we assumed in Eqns. (3.14) that the input vector is partitioned into 
the control input u z and the disturbance input wj. The software, however, can accommodate 
a more general form for the input vector. Using the variables ucol and ncol, defined later, 
the user can specify which inputs correspond to control and noise inputs respectively. Any 
row of the input vector can be either control, disturbance or both. (We interpret an input 
vector specified as both control and disturbance as two separate inputs, one a control input 
and the other a disturbance input whose distribution vectors have the same values.) 

Similarly, the outputs y c i and y are specified by crow and srow respectively, and 
any row of the output vector can be criterion output, sensor output or both. 

NOTE: For a discrete plant, the process noise distribution matrix is unity. 

wxx A vector of the PSD values of the white noise disturbance If 

Wfi(T) = E{wi(t- r)w[(t)\ (3.16a) 

orif W i (S(k)+8(j)) = E{w i (m,n)w'f (m + k,n + j)) (3.16b) 

depending on the value of stp , then 

wxx = [W-[ W2 ... Wjsjp] 


continued on following page... 
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Table 3.2 Matlab Workspace Variables Defining Optimization Problem (continued from previous page) 

Variable Description 

vr The PSD values of the discrete sensor noise covariance. If the output, sampled every STP, 

is 

>5#) = HiXi(k) + vi(k) (3.17) 

where V,8U) = E{vi(k)vf ( k + j)\ 

then 

YT=[V l V 2 ... V Np ] 


ncol A vector specifying which columns of G, correspond to If ncol = [1 4], then the plant 

disturbance vector w is comprised of the first and fourth plant inputs. 

ucol A vector specifying which columns of G( correspond to u r If ucol = [2 4], then the plant 

control input vector u is comprised of the second and fourth plant inputs. (Again, we 
interpret an input vector specified as both control and disturbance as two separate inputs, 
one a control input and the other a disturbance input whose distribution vectors have the 
same values.) 

srow A vector specifying which rows of Q correspond to y s j. If srow = [1 3], then the first and 

third continuous plant outputs are sampled and connected to the compensator. 

crow A vector specifying which rows of Q correspond to y Ci . If srow = [1 3], then the first and 

third continuous plant outputs are used to calculate the cost given in Eqn. (3.18). 


qla, q2a, maThe cost function weighting matrices for either a continuous or discrete cost function, 
depending on the value of stp defined later. The cost function has the form 


/ 


continuous 


-m 24s® 2K? 


(3.18a) 


or 


J discrete “ 


m — >00 


1 

Stppbtp , 


Np Stppbtp- 1 


1 = 1 n = 0 


jry CJ -(m,n)l 1 

r r0i« 

Mi' 



mJ 

Qii. 

Uf(m,n) JJ 


(3.18b) 


where i corresponds to a particular plant condition. 

qla =[Qn 012 - 0i Np] 
q2a = [<221 022 - 02Np] 
ma= [N] N 2 ... A^j] 

qa and ra must be symmetric and positive semi-definite. 

If ma is left undefined its value is assumed to be zero. 


continued on following page ... 
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Table 3.2 Matlab Workspace Variables Defining Optimization Problem (continued from previous page) 



cmp The compensator gain matrices as defined in Table 3.1. The digital processor gains can be 

either periodically time- varying or time -in variant. 

IMPORTANT: The user must provide an initial guess for cmp which stabilizes all plant 
perturbations. 

cmpfree This matrix specifies which of the gains in cmp are fixed, cmpfree has the same 

dimensions as cmp, but its elements are either 0 or 1 . AO indicates that the corresponding 
element in cmp is fixed and should not be optimized. A 1 indicates that the corresponding 
element in cmp is free and can be optimized. 

su, sy, sz The sampling schedule description. See Table 3.1. By selecting appropriate values for su, 
sy, sz and Stppbtp the designer can specify whether the compensator is multirate or single- 
rate. 

Stppbtp The same as stppbtp in Table 3,1. Note that the first letter of Stppbtp is upper case, 
indicating this is a global workspace variable. See Section 3.3.2. 

stp STP as defined in Table 3.1. 

If stp > 0 the algorithm assumes that the plant (pit), the processor noise (wxx) and the 
cost weighting matrices (qla, q2a, and ma) are all defined for a continuous 
plant in feedback with the discrete compensator. 

If stp = 0 the algorithm assumes that the plant (pit), the processor noise (wxx) and the 
cost weighting matrices (qla, q2a, and ma) are all defined for a discrete plant 
in feedback with the discrete compensator, and that the sampling period of the 
discrete plant is one STP of the multirate compensator 


3.3.2 mroptjnit 

Format : mropt_init 

Description : This script operates on the workspace variables defined in Table 3.2. It discretizes (if necessary) the 

plant, process noise and cost function; and defines a set of global workspace variables used by the 
optimization script (mropt_optim). This script must be rerun after making any changes in the 
workspace variables defined in Table 3.2, with exception to changes in the values of cmp or 

cmpfree. 

Inputs : mropt_init is a script which uses the workspace variables defined in Table 3.2. 

Outputs : mropt_init save it’ s output to the global workspace variables defined in Appendix A. 

Comment. The global workspace variables defined by mropt_init are used by both the optimization script 
mropt_optim and by the M-File calc_LQGcost . All global variables generated by mropt_init 
begin with a capital letter. This helps to differentiate them from other workspace variables. A 
brief description of these global variables is given in Appendix A. 
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3.3.3 mropt_optim 

Format : mropt_optim 

Description : mropt_optim computes optimum values of the digital processor gains such that the cost function 

defined in Eqn. (3.18) is minimized. 

Inputs : mropt_optim uses the global workspace variables generated by mropt_init in addition to the three 

workspace variables defined in Table 3.3. 

Outputs : The script returns a value of the cost function, the gradient of the cost function with respect to the 

free compensator parameters, and the optimum values of the digital processor gains. These 
values are automatically displayed on the computer screen and simultaneously stored in the 
workspace variables defined in Table 3.4. 

Comments'. 1) As already noted, if the user modifies the problem definition by changing the values of the 
variables defined in Table 3.2, the script mropt_init must be rerun before running mropt_optim . 
The exceptions are changes to the variables defined in Table 3.3. 

2) mropt_optini can be computationally intensive. The user can abort the optimization with a 
Control -C key sequence. Also see mropt_extract, Section 3.3.4 


Table 3.3 Input Workspace Variables for mropt_optim 



cmp This variable contains the compensator gains and is the same variable as described in 

Table 3.2. These are the starting values for the optimization and must stabilize all plant 
perturbations. 


cmpfree This matrix specifies which of the gains in cmp are fixed and is the same variable as 
described in Table 3.2. 

dscale During optimization the compensator gains are scaled using a variable named dscale to help 

improve the numerical accuracy of the search. They are scaled as follows 

Xscaled = mv(6scsdeyXunscaled 

where Xunscaled is an unsealed vector containing the free compensator gains and Xscaled 
are those values scaled. Typically dscale is a diagonal matrix chosen so that all the 
elements in Xscaled have the same magnitude. The user can either manually set dscale or 
have mropt_optim set it automatically. If dscale exists and has the correct dimensions 
then mropt_optim uses that value of dscale, otherwise a new value for dscale, based on the 
current values in cmp, is automatically calculated. Typing “clear dscale” before 
running mropt_optim forces the routine to scale the problem based on current values in 
cmp. 
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Table 3.4 Output Workspace Variables for mropt_optim 



cmpf This contains the optimized values of cmp after mropt_optim has completed the 

optimization. 

jcost This contains the value of the cost function after mropt_optim has completed the 

optimization. 

djdp This contains a scaled value of the gradient of the cost function with respect to the 

compensator values after mropt_optim has completed the optimization. 


3.3.4 mropt_extract 

Format: mropt_extract 

Description : After aborting the script mropt_optim, usually with a Control-C key sequence, the script 

mropt_extract can extract the value of the “optimized” compensator at the last iteration. 

Inputs : None. Can only be rerun after aborting mropt_optiin . 

Outputs : The optimized values of the digital processor gains at the last iteration before mropt_optim was 

aborted are displayed on the computer and stored in the workspace variable cmpf. 

Comments : Sometimes it is necessary to abort mropt_optim before is has completed the optimization. The user 

might wish to rescale the free parameters by changing dscale, or might decide to free up 
additional gains in cmp by changing cmpfree . The user can restart the optimization from the last 
iteration by setting cmp =cmpf and rerunning mropt_optim . 

3.3.5 calc_LQGcost 

Format : [jcost ]=calc_LQGcost(j) 

Description : Calculates the discrete LQG cost assuming a sampling period of one stp, This routine uses the 

global workspace variables computed by mroptjnit. Use this routine to find the minimum value 
of the cost function attainable by mropt_optim 

Inputs : j = 0 or is undefined, then the routine calculate the sum of the LQG costs for all Np plant conditions 

corresponding to the cost J discrete given in Eqn. (3.18). 

Ckj<Np then the routine calculates the LQG cost for only the plant condition 

Outputs : jcost The value of the LQG cost function 

Comments : calc_LQGcost uses the Controls_Toolbox routines dlqr and dlqe. They require that ra in Table 3.2 

be positive-definite. 

NOTE: Additional M-Files used in support of those described in this section are described in Appendix B. 
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SECTION FOUR 

Synthesis Example 


4.0 PROBLEM DESCRIPTION 

In this section we present an example design problem which utilizes the M-Files discussed in the previous 
section. Our objective is to design a multirate compensator for a lightly damped mass-spring-mass system. 

The mass-spring-mass (MKM) system is shown in Figure 4.1. It consists of two masses connected by a spring 
and damper. The control inputs are the force inputs u\ and U2\ the sensor outputs are the displacements x i and X2', 
and the disturbance input is w. The sensor outputs are corrupted by discrete sensor noise (not shown on figure) with 
covariance v. Nominal values for the plant parameters, along with those for one known plant perturbation, are given 
in Table 4.1. These parameter values result in a system with two poles at the origin associated with x\ , and two 
lightly damped poles associated with X 2 . The open-loop poles are listed in Table 4.1. A M-File which defines the 
system matrices for the MKM system is given in Appendix C. 

The goal of our design is to increase the damping coefficient on the lightly damped poles to £=0.707 without 
shifting their natural frequency, and to move the two poles at the origin so that the x\ response has a frequency of 
0.825 rad/sec with a damping coefficient £=0.707. 



Figure 4. 1 Mass- Spring-Mass System 


Table 4.1 Mass- Spring-Mass Parameter Values 



Nominal Plant 

Perturbed Plant 

M 

l.o kg 

1.0 kg 

m 

0.1 kg 

0.2 kg 

k 

0.01 N/m 

0.01 N/m 

b 

0.01 Ns/m 

0.01 Ns/m 

E{ w T (r-t)w(r)} 

0.16(t)N 2 

0.15(t)N 2 

E{v T (k+j)v(k)} 

0.001 S(/)m 2 

0.001 50')m 2 

Open-Loop Poles 

0, 0, -0.055+ 3.316/ 

0. 0, -0.030+2.449 i 
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4.1 COMPENSATOR DESIGN 

Our compensator design followed the preceding steps. All computer inputs and outputs in the following 
paragraphs are in this font. Bolded text should be input by the user, plain text is returned by the 
computer, instructional comments are in italic . 

Step 1: Select the weighting matrices for the quadratic cost function given in Eqn. (2.6) which characterize 

the desired closed-loop performance 

The cost function given in Eqn. (2.6) is the sum of the costs associated with each plant perturbation. In our case 
we have two plants - the nominal and the perturbed plant. We determined the weighting matrices for our cost 
function by designing two continuous time LQ regulators which placed the closed-loop poles for the nominal and 
perturbed plants in the desired locations. The cost functions are 

■^nominal = E{5.5x 2 +2.2x\ +10000u 2 + 10 k 2 }/10.6 (4.1a) 

/perturbed = e (6-5a: 1 2 + 4. 8i 2 + lOOOOu 2 + 10« 2 ) / 8. 0 (4.1b) 

where the costs have been scaled so that the LQR cost for each is unity. The cost function characterizing the 
performance of our multirate compensator is the sum of the these two cost functions. By minimizing (./nominal + 
♦/perturbed) optimization routine will attempt to find a single compensator whose performance approaches that of 
the LQR designs for both the nominal and the perturbed plant. Of course we do not expect such a compensator to 
exist. Instead the resulting compensator will represent a compromise between good performance at the nominal 
plant condition and good performance at the perturbed condition. 

Step 2: Select an appropriate compensator structure and sampling schedule based on the desired closed- 

loop performance 

The closed-loop LQR system has a fast mode associated with x\ and u\ and a slow mode associated with X2 
and W2* We selected a multirate compensator which capitalized on this structure. It consists of two coupled first- 
order loops, one from x\ to u\ and another X2 to «2 • A block diagram of this structure is shown in Figure 4.2. The 
x\ to wj loop is sampled and updated every 0.8 seconds while the X2 to «2 loop is sampled and updated four time as 
often, every 0.2 seconds. These sampling rates are approximately 10 time the desired closed-loop bandwith for each 
loop. In addition, the compensator accounts for a one-half sampling period computational delay. It’s sampling 
schedule is shown in Figure 4.3. 



x 1 Command 
*2 Command 


Figure 4.2 Block Diagram of Closed-Loop Multirate System 
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Y = sample/update activity 
Figure 4.3 Sampling Schedule for Multirate Compensator 


Step 3: Design a compensator with the sampling schedule and structure selected in Step 2 that stabilizes all 

plant perturbations. This will be the starting point for the optimization. 

We designed a compensator using successive loop closures and root locus for the nominal plant. The 
compensator has the same structure and sampling schedule as the multirate design discussed in Step 2 except there is 
no input coupling and we did not account for computational delay. This compensator does however stabilize both 
the nominal and perturbed plants even with the computational delay and so can be used as the starting point for the 
optimization. 

Refer to Reference 1 for an introduction to successive loop closures. 

Step 4: Load the problem definition into Matlab’s workspace variables defined in Table 3.2. 

We defined the workspace variables using the Matlab script provided in Appendix D. To load the workspace 
variables for this example, type 

»mropt_mtan 

Step 5: Initialize the workspace variables and generate the necessary global variable using mropt jnit 

Type 

»mropt_init see Section 3.3.2 

At this point you can compute the minimum value of the cost function as follows (see Section 3.3.5) 

»calc_LQGcost see Section 3.3.5 
» ans - 

2.0108e+00 this is the minimum cost our design could achieve 
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Step 6: Calculate optimum values of the digital processor gains using mropt_opdm 

A partial output of the optimization follows 

»mropt_optim see Section 3.3.3 

Automatically selecting dscale a diagonal dscale is automatically selected if the variable 

dscale does not exist, or is not compatible with the 
current compensator 

Calculate gradient only? (y or n)n if we replied y the algorithm would compute only 

the gradient of the cost with respect to the free 
digital processor gains and the value of the cost 
for the current compensator gains in cmp 

f -COUNT - the total number of function evaluations 
FUNCTION = the value of the cost function at the current iteration , Refer to Matlab FMINU 
documentation 

f -COUNT FUNCTION STEP-SIZE GRAD/SD LINE-SEARCH 

2 176 0.01 - 9 . 48e+07 

-> Plant 1 unstable while calculating J 

-> Plant 1 unstable while calculating J 

-> Plant 1 unstable while calculating J 

-> Plant 1 unstable while calculating J The algorithm encountered a destabilizing 

-> Plant 1 unstable while calculating J solution . It automatically adjusts the step 

-> Plant 1 unstable while calculating J size and continues the optimization 

-> Plant 1 unstable while calculating J 

-> Plant 1 unstable while calculating J 

-9 . 9655e-01 
-1 . 63 99e-02 
1 . 0030e+00 
-1 . 0676e-02 

-2 . 2977e-02 these are the values of the scaled digital processor gains 
-9.4363e-01 at the current iteration 
-1 . 3703e-03 
1 . 0529e+00 
1 . 003 0e+00 
1. 0036e+00 
1 . 0529e+00 
1 . 0085e+00 
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313 4.62 14.2 2.45e-13 int_st 

Optimization Terminated Successfully 
Gradient less than options (2) 

NO OF ITERATI0NS=3 14 
-8.0984e-01 
-1 . 67 00e+01 
1.8697e+00 
2 . 2849e+01 
7 . 4219e-01 
-7 . 3 479e + 00 
-8.6882e+00 
7 . 93 08e+00 
1 . 07 0 6e+00 
3 . 8424e~03 
1.0367 e+00 
1 . 4037e+00 


Final Results 


Gain 

Gradient 

-1.6197e-01 

-4 . 6527e-06 

-1.6700e+01 

6 . 6547e-04 

1.8697e+00 

-2 . 8675e-04 

2 . 2849e+01 

4 . 5 050e-04 

7 . 4219e-01 

1 . 657 6e-04 

-5 . 5109e+00 

-8 . 9955e-04 

-8 . 6882e+00 

~4 . 8353e-05 

7.930 8e + 00 

-8 . 5171e-04 

8 . 5648e-02 

-1 . 083 8e-04 

1 . 9212e-03 

6 . 9640e-04 

6 . 2202e-01 

3 . 4133e-03 

1 . 4037e-01 

3 . 2564e-04 


Final cost = 4.62482 the optimum value of the cost function 

Optimized compensator gains 
cmpf = 


1 . 9212e-03 

0 

1 . 8697e+00 

-8.6882e+00 

0 

1 . 4037e-01 

2 . 2849e+01 

7 . 9308e+00 

8 . 5648e-02 

0 

-1 . 6197e-01 

7.4219e-01 

0 

6 . 22 02e-01 

-1 .6700e+01 

-5.51 09e+00 


Elapsed time: 1516.05 sec 
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42 DESIGN ANALYSIS 

We looked at two criteria when evaluating our multirate design: 1) the final cost; and 2) the step response to a 
command input. 

4.2.1 Final Cost 

The final cost for the multirate compensator is calculated automatically by the optimization routine 
mropt_optim . The LQG cost is computed by the function calc_LQGcost in step 5. For our system 

^multirate — 4.6 and J LQR — 2.0 


4.2.2 Step Response 

We obtained the response to a unit command step to x\ as follows 
»[&,b,c,d] -mkm; get the nominal plant 

»b« [b zeros (4,l)];c*[c;[l 0 0 0]]; create a reference input 
»d*[d zeros (4,1);0 0 -1]; 

»ur* [ zeros (2,100); ones (1,100)]; generate a step input function for the input 

and compute the time domain response ; see Section 3.2.1 
» [y # t ] =mrsim( [a b;c d] , 4, cmpf , 2, su, sy, sz, Stppbtp, 1, [1 2] , [5 3],ur,stp); 

assuming cmpf su, sy, sz, Stppbtp, and stp were previously defined and 
that cmpf contains the optimized values of the digital processor gains 

»plot (t ,y ( [1 3],:)) 

We similarly obtained the step response for the LQR design. The results are shown in Figure 4.4. 

As expected the performance of the multirate design does not match that of the LQ regulator. The LQ regulator 
was optimized for only the nominal plant. The multirate designs on the other hand represent a compromise, 
stabilizing both plants and providing “optimum” performance for both. In addition, the LQ design was a continuous 
full-state design whereas the multirate design was discrete second order and used only two state measurements. 
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APPENDIX A 

Global Workspace Variables 


The following global workspace variables are defined by mropt_init Note that the first letter of every global 
variable is capitalized. 

For each plant perturbation, the discrete closed-loop system can be written in the form: 

Xi(m 1 ) = F(x{ m,n) + Gj u(m ,/i) + WW j (A. 1 ) 

yi(m,n)= Vw / (A. 2) 

Uj(m t n)= [S\(n)P(n)S2(n) + S3(n)]y s (m,ri) (A. 3) 


where 
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o' 


G l 

o‘ 


7 

o' 


Hi 

o" 


'0 

/ 


\D Z 

c z 

l 

, Gj = 

. w = 



, Hi- 



, v = 

0 

0 

, p = 

U 


_o 

0 

* i 

0 

/_ 

_0 

0_ 

’ l 

0 

/ 



A z. 


and w is a discrete zero mean white noise input with 

E{wj(m, n)wj ( m + j,n + fc)} = {<5(£) + 8(j)}Ri 


(A.4) 


(A.5) 


The subscript i indicates the ft plant condition and the overbarred matrices are the discretized plant matrices. The 
matrices 51, 52, and 53 are periodically time- varying switching matrices which model the multirate sampling 
scheme, and P is a matrix of the compensator’s digital processor gains. 

The discrete cost function for the system in Eqns. (A.1-A.3) can be written as 


J = lim 
m — ^°° 


± 

N 



Xi(m,n) 


iT 





(A.6) 


See Reference 1 for more information. 


The global Matlab workspace variables are then defined as follows. 


Plant Description: 

Np = number of plant conditions 


MF = [F\, F2, Fj sp] 
MG= [G], G2 , Gn p ] 
MW= [ W \ , W N p ] 

MV=[Vi,V 2 ,^ V Np ] 
MR= *Np] 


NcF 

NcG: 

NcW 

NcH= 

NcV= 

NcR = 


= number of columns in F( 

- number of columns in G( 
= number of columns in Wi 
number of columns in Hi 
number of columns in Vi 
number of columns in R( 
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Compensator Description: 

Nzhcmp = number of states in the periodically time- varying model of the multirate compensator 
Stppbtp = the global version of stppbtp 

Slk = [51(0), 51(0), 51(N-1)] Nsl = the number of columns in 51 

S2k = [52(0), 52(0), 52(N-1)] Ns2 = the number of columns in 52 

S3k = [53(0), 53(0), ..., 53(N-1)] Ns3 = the number of columns in 53 

Ptv: abs(Ptv) = number of columns in the compensator gain matrix [ A z , B z ; C Z> D Z ]. Ptv is positive if the 

digital processor gains are time-invariant and negative if they are periodically time-varying. 

Cost Function Description: 

MQ1= [fin, Q 12 , <2lNp] NcQl= number of columns in Qn 

MQ2= [ Q 21 , £>22* C?2Np] NcQ2= number of columns in 22,- 

MM= [M\ , M 2 , ... A^Np] NcM= number of columns in M/ 

Miscellaneous Definitions: 

Ppn , Ppm, and Ppinit : used internally 

Ppfree = index to the “free” digital processor gains 

Fd = closed-loop state transitions matrix of the last computed ETIS system 

Xc_global = value of the compensator gains at the last iteration 

Unstable : If Unstable =0 the compensator gains at the current iteration stabilize the plants 

Otherwise they destabilize the plant 
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Appendix B 

Support M-Files 


B.l MODELING SUPPORT M-FILES 


B.1.1 mrsplit 

Format : [a, b, c, d] = mrspiitQP, n) 

Description : Separates the composite gain matrix P into its constituents, where 


d c 

> = [b a 


(B.l) 


Inputs : P see Eqn. (B.l) 

n the number of columns in a 


Outputs : 


a,b,c,d the gain matrices in Eqn. (B.l) 


B.L2 mrmakesk 

Format : [slk, ns lk, s2k, ns2k, s3k, ns3k]= mrmakesk( su, sy, sz, stppbtp) 


Description : The periodically time-varying representation of a multirate compensator in Eqn. (2.2) can be written 

as 


y(m,/i) 
x(m,n + l) 


~D(n) C(n)Ju(m t n)l L.JW C z (n ) 

B(n) A(n)J[x(/n,n)J [ ^'[^(n) A^n) 


s: 


2(/i) + S3(n)j| 


u(m,n) 

x(m,n) 


(B.2) 


where A z , C z and D z are the digital processor gains and SI, 52 and 53 are switching matrices 
which model the compensator’s sampling schedule, mrmakesk creates the periodically time- 
varying switching matrices given the sampling schedule description. 

Inputs: su, sy, sz, stppbtp see Table 3.1 

Outputs: slk = [51(0) 51(1) ... Sl(stppbtp-l)] 

s2k = [52(0) 52(1) ... 52( stppbtp- 1)] 
s3k= [53(0) 53(1) ... 53( stppbtp- 1)] 

nslk, ns2k and ns3k the number of columns in 51(0, 52(0 and 53(0 respectively 


B.1J mrgetsk 

Format : [si , s2, s3]= mrgetsk (k, slk, nslk, s2k, ns2k, s3k, ns3k) 

Description: Extracts the individual switching matrices for a given STP from slk, s2k and s3k 

Inputs : k an integer specifying for which STP the switching matrices are to be extracted 

slk, nslk, s2k, ns2k, s3k, ns3k see Section B.l. 2. 

si, s2, s3 switching matrices corresponding to 51(k), 52(k) and 53(k) respectively. See 
Section B.l. 2. 


Outputs : 
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B.2 SYNTHESIS SUPPORT M-FILES 
B.2.1 calc_djdp 

Format : djdp=calc_djdp( pguess, dscale) 

Description : Used by the routine mropt_optim to calculate the derivative of the cost J with respect to the free 

compensator parameters. See Section 2.3. 

Inputs : pguess a vector of the values of the free digital processor gains at the current iteration 

dscale see Table 3.3. 

Comments: This routine uses the global workspace variables generated by mropt_init 

B.2. 2 calc_j 

Format : djdp= calcj (pguess, dscale) 

Description: Used by the routine mropt_optiin to calculate the cost J at the current iteration 

Inputs : pguess a vector of the values of the free digital processor gains at the current iteration 

dscale see Table 3.3. 

Comments: This routine uses the global workspace variables generated by mropt_init 

B.23 calc_L 

[lacc]= calc_L(pk, f, g, w, h, ql, q2, m) 

Called by calc_djdp to calculate the steady-state values of a Lagrange multiplier. (The multiplier is 
used to adjoin the stability equality constraints to the cost function.) 

pk gains for the periodically time- varying representation of the multirate compensator 

[D(i) C(i)l 

P k = [P(0),P(1),...,P(A^- 1)] and P(i) = (B.3) 

_B(i) A(i)_ 

D(i), C(i) t B(i), and A(i) are the state space matrices for the periodically time varying 
representation of the compensator. See Section 2.2 and Section B.1.2 
t g, w, h plant description matrices for the current plant perturbation corresponding to F Z f G;, W and 
H( in Eqn. (A.4). 

ql, q2, m cost weighting matrices for the current plant perturbation corresponding to Qn, £>2j and 
Mi in Eqn. (A.6). 

Output : lacc periodically time varying lagrange multipliers where 

lacc = [ A(l) A(2) ... A(A-1) A(0)] 
and A(0 is the lagrange multiplier for the i tfl STP 


Format: 

Description: 

Inputs : 


Comments: This routine uses the global workspace variables generated by mropt_init 
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B2A calc_P 

Format : [pk]= calc_P ( p) 

Description : Called by calc _j and calc_djdp to compute the periodically time-varying representation of the 

multirate compensator 

Inputs : p the digital processor gains in the form 

p = [p(0),p(l),...,p(A^-l)] where p(z) = 


D z (i) C z (i) 
B z {i) A z {i) 


Outputs : 


Comments: 


pk matrix of the compensator gains for the periodically time-varying representation of the multirate 
compensator 

pk = [P( 1), P(2\ P( 3),... P (N)} where F\i)= Sl(i)p(i)S2(i)+ S3(i) 

See Eqn. (A.2) or Section B.1.2. 

This routine uses the global workspace variables generated by mropt _init 


B.2.5 calc_X 


Format : 
Description : 

Inputs : 


Comments : 


[xacc]- calc_X(pk, f, g, w, h, v, r) 

Called by calc_djdp and calcj to calculate the steady-state covariance values of the closed-loop 
system 

pk matrix of the compensator gains for the periodically time-varying representation of the multirate 
compensator, see Section B.2.4. 

£, g, w, h plant description matrices for the current plant perturbation corresponding to > G/, W, H( 
and V in Eqn. (A.4) respectively 

r process and sampling noise covariances for the current plant perturbation corresponding to R( in 
Eqn. (A.5). 

This routine uses the global workspace variables generated by mropt_init 


B.2.6 disc_cost 


Format: 

Description: 


[qld, q2d , md]=disc_cost (pit, nplt, stp, ql, q2, m) 


Computes the discrete cost function weighting matrices such that the cost associated with a 
continuous plant in feedback with a single-rate compensator is identical to the cost associated 
with the discretized plant in feedback with the single-rate compensator. The continuous cost 
function is given in Eqn. (2.6); the discrete cost function is given by 


J discrete “ \ j 


N-l 


x* 


fr V 

| x(m,n) 

r r Qid 

M d ' 

rx(m,n) 

u(m,«) 

mJ 

Qid. 

\u(m,n) 


(B.4) 


disc_cost assumes the control input is updated using a zero-order-hold. 

Inputs: pit, nplt, stp see Table 3.1 

ql, q2, m the continuous cost function weighting matrices in Eqn. (2.6) 

Outputs: qld, q2d, md the discrete weighting matrices in Eqn. (B.4) 
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B.2.7 disc_noise 

Format : [w}= disc_noise(wxx, pit , nplt, stp ) 

Description: Computes the covariance of a discrete-time white noise process disturbance such that the expected 

value of the states of a continuous plant, and the expected value of the states of the discretized 
plant are equivalent at the sampling instances. Thus, given the continuous plant 

i(f) = Ar(f) + G w w(t), where E(w(t)w T (t + t)) = <5(t)wxx 

and the plant discretized with a zero-order-hold at a sampling period T 

xj(m y n + 1) = A^x^m.n)^ H^(/w,n), where E(wy(m,n)wj (m + k y n + j)) = (S(k) + <5(j))w 

then E{x^(m } n)xJ(m y n)} = E{x((mN + n)T)x^ ((mN + n)T)} 

Inputs : wxx PSD of the continuous processor noise 

pit, nplt, stp see Table 3.1. The inputs to pit are assumed to be only the process noise 

Outputs : w PSD of the discrete process noise 

B.2.8 dlyap2 

Format : x=dlyap2(a, c, method ) 

Description : Solves the discrete lyapunov equation 

x=a-x*a T +c (B.5) 

Inputs : a, c system matrices in Eqn. (B.5) 

if method = ‘Bartels’ then the algorithm solves Eqn. (B.5) using Bartels method 
otherwise Eqn. (B.5) is solved using partial sums 

Outputs : x the solution to the lyapunov equation 

B.2.9 get_cost 

Format : [ql , q2, m]=get_cost(i) 

Description : Extracts the cost function weighting matrices Ql, Q2 and M for the i th plant perturbation from the 

global workspace variables MQ1 , MQ2, and MM. See Appendix A. 

Inputs : i an index to the desired plant perturbation 

Outputs : ql , q2, m the cost weighting matrices corresponding to Ql, Q 2 and M respectively 

This routine uses the global workspace variables generated by mropt_init 


Comments : 
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B.2.10 get_plt 


Format : 
Description : 


[f, g, w, h, v, r)=get_plt(i) 

tUo. V. n ■ TJ7 XJ. T/ n-ryA D ■ fr\ r tV>o l th nlont n^rtirrViQti on frnm a1r»Ka1 

xtracts me plant matrices i j, Uj, r > , , v , anw * vj iui 

workspace variables MF, MG, MW, MH, MV and MR. See Appendix A. 


Inputs : i an index to the desired plant perturbation 

Outputs : f, g, w, h, v, r plant description matrices for the i^ 2 plant perturbation corresponding to F,\ G[, W Hi » 

V, and /?| respectively 


Comments : This routine uses the global workspace variables generated by mropt.init 


B.2.11 get_jppfree 

Format : ppfree = gefrppfree (empfree , nemp ) 

Description: Called by mropt.optim to determine which compensator gains will be optimized 

Inputs : empfree , nemp see T able 3.1. 

Outputs : ppfree vector whose elements point to the free elements in emp 


B.2.12 get.sk 

Format : [si , s2 , s3 ]= get_sk( k) 

Description: Extracts the switching matrices Sl(k ), 52(k) and S3(k) for the k^ 1 STP from the global workspace 

variables Slk, S2k, and S3k. See Appendix A. 

Inputs : k indicates the k^ 1 STP of the current BTP. NOTE: k = 1, 2, ... , stppbtp 

Outputs : si, s2, s3 switching matrices for the k ^ STP corresponding to 51 (k), 52(k) and 53(k) respectively 

Comments : This routine uses the global workspace variables generated by mropt.init 

B.2.13 make.noise 

Format : [r ]= make.noise (rw , rv) 

Description : Creates a compound noise covariance matrix from rw and rv 

Inputs : rw , rv the process and sensor noise covariance matrices respectively 

Outputs : r compound noise matrix of the form 

rw 0 
0 rv 
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B.2.14 make^plt 

Format : [f, g, w, h, v, nzcmp >= make_plt (pltc, nplt, ncmp, stp) 

Description : Creates the compound plant matrices Fj, G* } R 7 and 7 for the current plant perturbation. See 

Appendix A. 

Inputs : pltc has the same form as pit in Table 3.1 except it describes only a single plant perturbation 

nplt, ncmp , stp See Table 3.1. 

Outputs : f, g, w , h, v compound plant matrices F*, G,\ W and V 7 for the current plant perturbation respectively 

nzcmp the number states in the periodically time-varying representation of the multirate 
compensator 

B.2.15 mropt_check 
Format: mropt_check 

Description: mropt_check performs elementary error checking on the data in Table 3.2 and set the variable 

err=l if it finds an error in the data. mropt_check is called by mropt_init 

Inputs : None 

Outputs : None 

B.2.16 mropt_global 

Format : mropt ^global 

Description: Called by mropt_init to define the global workspace variables 

Inputs: None 

Outputs : None 

B.2.17 mroptfminu 

Format: [x, OPTIONS] = mropt_fminu(FUN, x, OPTIONS, GRADFUN, PI, P2, P3, P4, P5, ... 

P6, P7, P8, P9, P10) 

Description : A modified version of FMINU from the Optimization_.Toolbox. This version automatically reduces 

the search step size when a destabilizing solution is encountered. 

Inputs : See FMINU 

Outputs : See FMINU 

Comments: Requires the Optimization_Toolbox 


i 
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APPENDIX C 


M-Filemkm 


The following M-File creates the state space matrices for the Mass-Spring-Mass system of Section 4.0 

function [a, b, c f d] =mkm (ml .m2 , k, b) 

% [a, b, c, d] =mkm{ml , m2 , k, b) 

if nargin~=4 %nominal plant 
ml=l ; 
m2= . 1; 
k=l; 
b=0 . 01 ; 
end; 

t= (l/ml)+{l/m2) ; 

a= [ 0 1 0 0 

0 0 k/ml b/ml 

0 0 0 1 
0 0 -k*t -b*t] ; 

b= [0 0 

1/ml 0 
0 0 

-1/ml l/m2]; 

c=eye (4) ; 
d=zeros (4,2) ; 
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Blank 
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APPENDIX D 

Script mropt_mkm 

The following script defines the workspace variables in Table 3.2 for the example problem of Section 4.0 
% Script: mropt^mkm . m 

% Example script for creating input data for mropt 
% 

% Creates the following data required by mropt_init 
% plt,nplt ,wxx, rv, qla, q2a,ma 

% ucol , ncol , srow, crow 
% sz f su, sy, cmp, stp, Stppbtp, cmpfree 

% G. Mason, U.W. 1992 

% Continuous plant description 

Np=2; % The number of plant conditions 

nplt=4; % The number of states in the plants 

plt=[]; % Clear the variable which holds the plant data 

% Plant #1 

[a, b, c, d] =mkm(l , 0.1,1,0.01) ; 
pit = [a b; c d] ; 

% Plant #2 

[a,b, c , d] =mkm{ 1 ,0.2,1,0.01); 
pit = [pit [a b ; c d] ] ; 

% Pointer to rows and columns in the plant 
ucol=[l 2]; % control input 

ncol=2; % process noise input 

srow=[l 3]; % sensor output from the plant 

crow=[l 234]; % criterion output from the plant 

% Continuous process noise PSD 

% wxx is a coirpound matrix like pit, it contains wxx for plant #1 and #2 
wxx= [ . 1 . 1 ] ; 

% Discrete sensor noise PSD 
% with same compound form as pit 
rv=eye (2 ) *0 . 001 ; 
rv= [rv rv] ; 

% Continuous Cost weighting matrices 

% Output criterion weights 

% Use le-8 instead of 0 so that qal & qa2 are positive definite 
% The synthesis algorithm would accept semi-definite qla & q2a 
% but calc_LQGcost will not 

qlal=diag( [5.5 le-8 le-8 2.2]); 
qla2=diag ( [ 6 . 5 le-8 le-8 4.8]); 
qla= [qal/10 . 6 qa2/8 . 0] ; 

% Control input weights 
q2a=diag( [le4 lei] ) ; 
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q2a= [ra/10 . 6 ra/8 . 0] ; 

% Cross weighting (ma) 

% if there is no cross weighting it can be left undefined 
clear ma 


% Compensator description 

% Digital processor gains from successive loop closures design 
cmp= [0.50 1 0 

0 0.1 0 1 
0.08 0 -0.2 0 
0 0.6 0 -0.75]; 


% Sampling schedule 

su= [ 0 0 ; 1 1 ; 0 0;0 1;0 0;0 1;0 0;0 1] ; 

% Compensator output updating w/ delay 
sy=[8 2]; % Sensor input sampling, multiplexed 

sz=[8 2]; % Compensator state update 

stp = .1; % The shortest sampling period 

Stppbtp =8; % The number of stp’s in one BTP 

% Free compensator gains 
cmpfree=[l Oil 
0 111 
10 11 
0 111]? 


% clear temporary variables 
clear abed 
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